#!/usr/bin/env Rscript
#######################
## Interaction plots
## For fraud project
#######################

# Flag to send job to cluster or 
# run in a regular computer (still in parallel,
# using MPI backend)
PegasusRun  <- commandArgs()

if(PegasusRun){
  system("bsub < BartIneractionPlots.job")
}else{
  
  ## Load required libraries
  library(plyr)
  library(doMPI)
  library(BayesTree)
  
  ## Seed RNG
  set.seed(831213)
  
  
  ## Start cluster and worker loop
  cl  <- startMPIcluster()
  registerDoMPI(cl)
  
  
  ## Get processed data
  load(".-FullData.RData")
  train.data  <- all.data$fData[,-grep("uid|train|fraud",colnames(all.data$fData))]
  train.y  <- all.data$fData[,"fraud.score"]
  
  
  ## Automatically obtain levels at which
  ## partial dependence plot is evaluated.
  ## Variable indeces are based on column numbers 
  ## of train.data as defined on line 23 above.
  inter.arguments <- list(
    list(xind=c(7,2)
         ,levs=list(seq(min(train.data[,"turnout"]),max(train.data[,"turnout"]),length.out=11)
                    ,seq(min(train.data[,"DistanceLastPair"]),max(train.data[,"DistanceLastPair"]),length.out=11))
    )
  )
  
  ## Get best parameters from parameter sweep
  bestParams <- read.csv("./BSParamSweep.csv"
                         ,nrows=1
  )[,1:6]
  
  ## Use pd2bart to estimate two-way partial dependencies
  bartInterPlots <- llply(inter.arguments
                          ,.fun=function(args){
                            pd2bart(x.train = train.data
                                    ,y.train = train.y
                                    ,xind=args[[1]]
                                    ,levs=args[[2]]
                                    ,pl=FALSE
                                    ,k = bestParams$k
                                    ,power = bestParams$pow
                                    ,base = bestParams$base
                                    ,sigdf = bestParams$sigdf
                                    ,sigquant = bestParams$sigquant
                                    ,nskip = 50e3
                                    ,ndpost = 5e3
                                    ,keepevery=10
                                    ,keeptrainfits=FALSE
                                    ,ntree = 100
                                     ,verbose = FALSE
                            )                          
                          }
                          ,.parallel=TRUE
                          ,.paropts=list(.packages="BayesTree"
                                         ,.export=c("train.data","train.y"))
  )
  
  save(bartInterPlots,file="./BartPlotObjectsInter.RData")
  
  closeCluster(cl)
  mpi.finalize()
}
